--- title: "Class 3" author: "Colin Kuehl" date: "`r format(Sys.time(), '%B %d, %Y')`" output: html_document editor_options: chunk_output_type: console --- This week we'll start by illustrating some probability concepts in R and then continue building R skills by looking at the gapminder dataset # Part 1: Probability ######################################################### ``` ## R code to demonstrate probability distributions. You won't need to use this code later so don't worry if you are unclear what is going on. ``` Brief illustration of the law the large numbers ```{r} #Don't worry about the code here. This is just an illustration. n<-c(1:100) #creates a vector of 1 to 100 weight<-rnorm(100,3,.5) #Create a vector of length 100, with a mean of 3, and sd of .5 hist(weight)#quick histogram of what we just created s<-cumsum(weight) #cumulative sum wmean <- s/n #finds mean at each point plot(wmean,xlab="n",ylim=c(2,4),type="l",main="Law of Large Numbers") abline(h=3, col="red")#horizontal line at the mean #what happens if we change to a thousand ``` ## simulate die rolls ```{r} d1 <- sample(1:6, 100000, prob=rep(1/6, 6), replace=TRUE) d2 <- sample(1:6, 100000, prob=rep(1/6, 6), replace=TRUE) hist(d1)#Expected uniform distribution for a singe die die.roll <- d1 + d2 #sum of 2 dice hist(die.roll, breaks= seq(0, 12, by=.5), prob=FALSE, col="red", main="Frequency Dist, 100,000 die rolls", xlab = "Value of Roll", ylab="Frequency") summary(die.roll) table(die.roll) prop.table(table(die.roll)) ## plot PMF freq<-cbind(table(die.roll)) prval<-freq/100000 roll<-cbind(2,3,4,5,6,7,8,9,10,11,12) plot(roll, prval, col="red", main="Probability Mass Function", xlab = "Value of Roll", ylab="Probability", pch=19, type="l") ## plot CDF #CDF = cumulative distribution function - probability that X will take a value less than or equal to X cprval <- cumsum(prval) plot(roll, cprval, col="red", main="Cumulative Mass Function", xlab = "Value of Roll", ylab="Probability", pch=19) ``` Creating different distributions ```{r} ## The uniform distribution x <- runif(10000) plot(x, x/x, type="l", col="red", main="PDF of Uniform", xlab = "X", ylab="Probability") lb <- 0.25 ub <- 0.75 polygon(c(lb,lb,ub,ub), c(0,1,1,0), col="red") ## The normal distribution x <- sort(rnorm(10000)) plot(x, dnorm(x), type="l", col="red", main="PDF and CDF of Normal", xlab = "X", ylab="Probability", ylim=c(0,1)) lines(x, pnorm(x)) abline(v=0) ``` ################################################################################################## #Packages This week we'll use some packages. Packages are mini pieces of software(sets of commands) that are loaded on top of base R. There are thousands to choose from and are very helpful for doing more specialized types of analysis. I recommend loading all the packages you need at the beginning of your code (along with setting your working directory and loading your data). The install command downloads the package to your computer. Once installed you don't need to reinstall every time. (in fact, sometimes its better not to in case of differences between versions) Two of the most commonly used are dplyr and ggplot The library command unpacks(makes available) the package. Like taking it out off the shelf. ```{r} #install.packages("tidyverse") ## This package is used to tidydata. We'll use it later #install.packages("ggplot2") #Very popular graphing package library(tidyverse) library(ggplot2) #I use these two packages almost every R session #setwd("/Users/colinkuehl/Dropbox/POLS 641/F26 641/f26course/Rcode/Week 3") ``` Today we'll explore the gapminder data. It contains life expectancy, population, and GDP per capita for countries every five years from 1952 to 2007. The gapminder data comes from the Gapminder Foundation: made famous by Hans Rosling's TED talks. We'll continue our exploration using R to understand and visualize data. ## Part 1: Getting to Know Gapminder Make sure you have opened the Rproject first. Today's data comes packaged inside R itself. Instead of loading a file, we install a package once, then load it with library() each session. Packages are like bonus material for R. We will use them all the time. Installing puts it on your computer. Library opens it(like a book) so everytime you use it you need library in your code first. ```{r} #install.packages("gapminder") #only run this once ever - then put the # back in front library(gapminder) data("gapminder") #Puts into our environment ``` Looking over the data using commmands from previous weeks (hint: use code from yesterday) ```{r} head(gapminder) tail(gapminder) summary(gapminder) dim(gapminder) #Count of how many observations and variables str(gapminder) #This tells us what kind of data R thinks each variable is - we'll talk about this later length(unique(gapminder$country)) #gives a count of unique countrys ``` Checking in (with a partner) 1) How many countries are in the dataset? 2) What is the average life expectancy? 3) What is the median GDP per capita? Average? Bonus: a) Create a scatterplot with the relationship between population and Life Expectancy. Add a trend line b)What country has the highest life expectancy? ```{r} plot(gapminder$lifeExp, gapminder$pop) table(gapminder$year) table(gapminder$country) ``` *It is best practice to "knit as you go." Try knitting now to make sure you don't have errors in your code.* Country and continent are not continuous(numeric in R speak) variables. They are categorical. We can get counts in each category. frequency tables: ```{r} table(gapminder$continent) #gives count of different categories table(gapminder$year) prop.table(table(gapminder$continent)) #gives percentages ``` I use this command all the time to get a sense of the data. Check in: What percentage of the observations in the dataset are from Asia? ```{r} ``` `` Plot the relationship between GDP per capita and life expectancy ```{r} plot(gapminder$gdpPercap, gapminder$lifeExp) #scatterplot (X then Y) abline(lm(gapminder$lifeExp~gapminder$gdpPercap)) #add trend line (Y then X) ``` Notice how we always put the explanatory variable on the X axis(horizantal) and the outcome variable on the Y axis (vertical) GDP per capita hard to read - change unit to thousands ```{r} gapminder$gdpPercap <- gapminder$gdpPercap/1000 ``` Make plot look a bit more professional ```{r} plot(gapminder$gdpPercap, gapminder$lifeExp, main="Relationship between Wealth and Life Expectancy", xlab="GDP per Capita", ylab="Life Expectancy") abline(lm(gapminder$lifeExp~gapminder$gdpPercap)) #and boom, you've already ran your 2nd regression ``` Does there appear to be a relationship in the sample? What about with just Asia? There are multiple ways to do this: We can create a new dataset ie a subset(I almost always choose this option) ```{r} table(gapminder$continent) asia <- gapminder[which(gapminder$continent=="Asia"),] asia <- filter(gapminder, continent=="Asia") # This does the same thing using the dplyr package from tidyverse. I find it much more intuitive than "base R" so tend to use it more often. plot(asia$gdpPercap, asia$lifeExp, main="Relationship between Wealth and Life Expectancy (Asia)", xlab="GDP per Capita", ylab="Life Expectancy") abline(lm(asia$lifeExp~asia$gdpPercap)) #and boom, another regression ``` Exclude Europe ```{r} noteuro <- gapminder[which(gapminder$continent!="Europe"),] # the != does not equal noteuro <- filter(gapminder, continent!="Europe") ``` Exclude the super rich - use greater than or less than ```{r} notrich <- gapminder[which(gapminder$gdpPercap< 50),] # less than 50,000 - this drops oil-rich outliers like Kuwait notrich <- filter(gapminder, gdpPercap< 50) ``` Check in: Try creating a scatterplot with just the 1) Africa OR 2) country-years where life expectancy is less than 70. ```{r} ``` #Looking at the data over time: Time to look at something a bit more fun. How has life expectancy changed across the world since the 1950s? Mean life expectancy? ```{r} mean(gapminder$lifeExp) min(gapminder$lifeExp) max(gapminder$lifeExp) range(gapminder$lifeExp) ``` Create a histogram for life expectancy ```{r} hist(gapminder$lifeExp) abline(v=mean(gapminder$lifeExp), col="blue", lwd=2)#adds a vertical(v) line for the mean abline(v=median(gapminder$lifeExp), col="red", lwd=2)#a a line for mediant ``` abline gives us a line. v is for vertical Density plots and boxplots show us similar information ```{r} boxplot(gapminder$lifeExp) plot(density(gapminder$lifeExp)) #How can we find out which country-years are at the very bottom? ``` Check: How many observations are from Africa and how many are from Europe ```{r} ``` ```{r} mean(gapminder$lifeExp[gapminder$continent=="Africa"]) #average for Africa mean(gapminder$lifeExp[gapminder$continent=="Europe"]) #average for Europe mean(asia$lifeExp) #find the mean for the Americas ``` Let's make a simple scatterplot with year on the horizontal axis, and life expectancy on the vertical axis. A regression allows us to look at the relationship between variables linear regression model, and you're looking at the effect of the IV on the DV ```{r} plot(gapminder$year, gapminder$lifeExp, ylab = "Life Expectancy", xlab = "Year", main = "Life Expectancy Over Time") abline(lm(gapminder$lifeExp ~ gapminder$year), col="red", lwd=2) # What type of relationship is this? abline(h=mean(gapminder$lifeExp), col="blue", lwd=2) ``` That is a bit overwhelming lets just look at 2007 - ie lets create a subset ```{r} g2007 <- gapminder[which(gapminder$year==2007),] #just using base R #New Graph plot(g2007$gdpPercap, g2007$lifeExp, ylab = "Life Expectancy", xlab = "GDP per Capita", main = "Wealth and Life Expectancy in 2007", col=factor(g2007$continent)) #color for continent abline(lm(g2007$lifeExp ~ g2007$gdpPercap), col="red", lwd=2) legend("bottomright", legend = levels(factor(g2007$continent)), pch=19, col = factor(levels(factor(g2007$continent)))) #I want to add labels for specific countries g2007$stars <- ifelse(g2007$country=="United States"|g2007$country=="China" |g2007$country=="Myanmar" |g2007$country=="Thailand", as.character(g2007$country), "") #create a new variable for notable countries, all else will be empty. text(g2007$gdpPercap, g2007$lifeExp-1, g2007$stars, cex=.4) #add text to location ``` Check: Can you create the same for 1952 ```{r} ``` A boxplot for each year can show us how the distribution has changed over time ```{r} boxplot(gapminder$lifeExp~gapminder$year, xlab = "Year", ylab="Life Expectancy", border = "steelblue") boxplot(gapminder$lifeExp~gapminder$continent, xlab = "Year", ylab="Life Expectancy", border = "steelblue") ``` Lets look at changes in thailand over time ```{r} thai <- filter(gapminder, country=="Thailand") plot(thai$year, thai$lifeExp, xlab="Year", ylab="Life Expectancy", col="red", main="Thailand Life Expectancy") lines(thai$year, thai$lifeExp, col="blue") ``` ```{r} plot(asia$year, asia$lifeExp, xlab="Year", ylab="Life Expectancy", col=gapminder$country, main="Thailand Life Expectancy") lines(thai$year, thai$lifeExp, col="blue") ``` Check: Do the same above for a different country I've clearly had too much fun. Your turn to make something interesting. ```{r} ```